/*
 * Multiphase.cpp
 *
 *  Created on: 23 Jul 2011
 *      Author: Allan
 */

#include "Multiphase.h"

// C++ includes
#include <algorithm>
#include <iomanip>
#include <tuple>

Multiphase::Multiphase()
{}

Multiphase::Multiphase(const vector<Phase>& phases) :
phases(phases)
{
	// Loop over all the phases to set the sequence of species 
	BOOST_FOREACH(const Phase& phase, phases)
		species.insert(species.end(), phase.GetSpecies().begin(), phase.GetSpecies().end());
	
	speciesMap = SpeciesMap(species);
}

const vector<Phase>& Multiphase::GetPhases() const
{
	return phases;
}

const vector<string>& Multiphase::GetSpecies() const
{
	return species;
}

const SpeciesMap& Multiphase::GetSpeciesMap() const
{
	return speciesMap;
}

const Index Multiphase::operator[](const string& species) const
{
	return speciesMap[species];
}

const vector<Index> Multiphase::operator[](const vector<string>& species) const
{
	return speciesMap[species];
}

template<class PhaseFunction> 
const VectorXd GenericPhaseComputation(double T, double P, const VectorXd& n, const vector<Phase>& phases, const PhaseFunction& phaseFunction)
{
	// The activity-like data of the species to be computed
	VectorXd a(n.rows());
	
	// The current number of processed species, to be update at the end of every loop (see below)
	unsigned nprocessed = 0;
	
	// Loop over all the phases
	BOOST_FOREACH(const Phase& phase, phases)
	{
		// The number of species in the current phase
		const unsigned nspecies = phase.GetSpecies().size();
		
		// The vector of number of moles of the species of the current phase
		const VectorXd np = n.segment(nprocessed, nspecies);
		
		// Compute the activity-like data of the species in the current phase using the template parameter phaseFunction
		const VectorXd ap = phaseFunction(phase, T, P, np);
		
		// Transfer the activity-like data of the current phase to the vector of activity-like data a (beginning at position given by nprocessed)
		a.segment(nprocessed, nspecies) = ap;
		
		// Update the number of processed species
		nprocessed += nspecies;
	}

	return a;
}

const VectorXd Multiphase::ActivityCoefficients(double T, double P, const VectorXd& n) const
{
	auto activitycoeffsFunc = 
			[](const Phase& phase, double T, double P, const VectorXd& n) { return phase.ActivityCoefficients(T, P, n); };
	
	return GenericPhaseComputation(T, P, n, phases, activitycoeffsFunc);
}
	
const VectorXd Multiphase::Activities(double T, double P, const VectorXd& n) const
{
	auto activitiesFunc = 
			[](const Phase& phase, double T, double P, const VectorXd& n) { return phase.Activities(T, P, n); };
	
	return GenericPhaseComputation(T, P, n, phases, activitiesFunc);
}

ostream& operator<<(ostream& out, const Multiphase& multiphase)
{
	out << "Multiphase:" << endl;
	BOOST_FOREACH(const Phase& phase, multiphase.GetPhases())
	{
		out << '\t' << "Phase: " << phase.GetName() << endl;
		
		BOOST_FOREACH(const string& s, phase.GetSpecies())
			out << "\t\t" << "(" << multiphase[s] << ") " << s << endl;
		
		out << endl;
	}
	
	return out;
}
